﻿using System;
using System.Collections.Generic;
using System.Text;

namespace MLForgeSharp.Models.ProbabilisticModels
{
    /// <summary>
    /// 高斯混合模型
    /// </summary>
    public class GaussMixModel
    {
        private int numComponents; // 高斯分布的数量
        private double[] weights; // 混合系数
        private double[] means; // 均值
        private double[] variances; // 方差

        public GaussMixModel(int numComponents)
        {
            this.numComponents = numComponents;
            this.weights = new double[numComponents];
            this.means = new double[numComponents];
            this.variances = new double[numComponents];

            // 初始化参数
            Random rand = new Random();
            for (int i = 0; i < numComponents; i++)
            {
                weights[i] = 1.0 / numComponents; // 均匀分布
                means[i] = rand.NextDouble() * 10; // 随机均值
                variances[i] = rand.NextDouble() * 10; // 随机方差
            }
        }

        // 计算高斯分布的概率密度函数
        private double GaussianPDF(double x, double mean, double variance)
        {
            return Math.Exp(-0.5 * Math.Pow((x - mean), 2) / variance) / Math.Sqrt(2 * Math.PI * variance);
        }

        // E 步：计算每个数据点属于每个高斯分布的概率
        private double[,] ExpectationStep(double[] data)
        {
            int numDataPoints = data.Length;
            double[,] responsibilities = new double[numDataPoints, numComponents];

            for (int i = 0; i < numDataPoints; i++)
            {
                double sum = 0.0;
                for (int j = 0; j < numComponents; j++)
                {
                    responsibilities[i, j] = weights[j] * GaussianPDF(data[i], means[j], variances[j]);
                    sum += responsibilities[i, j];
                }
                for (int j = 0; j < numComponents; j++)
                {
                    responsibilities[i, j] /= sum; // 归一化
                }
            }

            return responsibilities;
        }

        // M 步：更新模型参数
        private void MaximizationStep(double[] data, double[,] responsibilities)
        {
            int numDataPoints = data.Length;

            for (int j = 0; j < numComponents; j++)
            {
                double sumResponsibilities = 0.0;
                double sumMean = 0.0;
                double sumVariance = 0.0;

                for (int i = 0; i < numDataPoints; i++)
                {
                    sumResponsibilities += responsibilities[i, j];
                    sumMean += responsibilities[i, j] * data[i];
                }
                means[j] = sumMean / sumResponsibilities;

                for (int i = 0; i < numDataPoints; i++)
                {
                    sumVariance += responsibilities[i, j] * Math.Pow(data[i] - means[j], 2);
                }
                variances[j] = sumVariance / sumResponsibilities;

                weights[j] = sumResponsibilities / numDataPoints;
            }
        }

        // EM 算法：训练模型
        public void Train(double[] data, int maxIterations = 100)
        {
            for (int iter = 0; iter < maxIterations; iter++)
            {
                // E 步
                double[,] responsibilities = ExpectationStep(data);

                // M 步
                MaximizationStep(data, responsibilities);
            }
        }

        // 模型参数
        public string GetParameters()
        {
            string result = "";
            for (int j = 0; j < numComponents; j++)
            {
                result += string.Format("Component {0}: Weight = {1:F4}, Mean = {2:F4}, Variance = {3:F4}\n",
                                        j + 1, weights[j], means[j], variances[j]);
            }
            return result;
        }
    }

    public class GaussMixModelExample
    {
        public GaussMixModelExample()
        {
            // 示例数据（从两个高斯分布生成）
            double[] data = new double[] { 1.2, 1.8, 2.3, 3.5, 4.1, 5.0, 5.5, 6.0, 6.5, 7.0 };

            // 创建高斯混合模型
            GaussMixModel gmm = new GaussMixModel(numComponents: 2);

            // 训练模型
            gmm.Train(data);

            // 打印模型参数
            gmm.GetParameters();
        }
    }
}
